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ABSTRACT 

Analytical approaches to galaxy formation and reionization are based on the mathe- 
matical problem of random walks with barriers. The statistics of a single random walk 
can be used to calculate one-point distributions ranging from the mass function of viri- 
alized halos to the distribution of ionized bubble sizes during reionization. However, 
an analytical calculation of two-point correlation functions or of spatially-dependent 
feedback processes requires the joint statistics of random walks at two different points. 
An accurate analytical expression for the statistics of two correlated random walks has 
been previously found only for the case of a constant barrier height. However, calculat- 
ing bubble sizes or accurate statistics for halo formation involves more general barriers 
that can often be approximated as linear barriers. We generalize the two-point solution 
with constant barriers to linear barriers, and apply it as an illustration to calculate 
the correlation function of cosmological 21-cm fluctuations during reionization. 

Key words: galaxies:high-redshift - cosmology:theory - galaxies:formation - large- 
scale structure of universe - methods: analytical 



1 INTRODUCTION 

A critical prediction of any theory of structure formation 
is the mass function of virialized dark-matter halos. While 
only numerical simulations capture the full details of halo 
collapse, much of our understanding of structure formation 
relies instead on analytical techniques. As such methods 
are based on simple assumptions and are easily applied to 
a large range of models, they are indispensable both for 
gaining physical understanding into the numerical results 
and exploring the effects of model uncertainties. Analyti- 
cal methods also can be used to study and compensate for 
various limitations of numerical simulations such as insuffi- 
cient small-scale resolution, missing large-scale fluctuations, 
or insufficiently early starting redshifts. 

The mos t widely applied m e thod of this type was first 
developed by iPress fc Schechterl (Il974h. This simple mode l, 
later refined by iBond et al.1 (|l99ll ). lLacev fc Cole] (|l993h . 
and others, has had great success in describing the forma- 
tion of structure, reproducing rather accurately the numer- 
ical results. Yet this model is intrinsically limited since it 
can only predict the average number density of halos. Bary- 
onic objects forming within these halos are often subject to 
strong environmental effects that are untreatable in this con- 
text. Many environmental effects such as photoionization or 
metal enrichment are highly inhomogeneous in nature, be- 
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ing caused by the nonlinear structures that form within the 
intergalactic medium (IGM), and thus primarily impacting 
the areas near these structures. Such interactions between 
the IGM and structure formation are often better described 
as spatially-dependent feedback loops rather than sudden 
changes in the overall average conditions. 

The issue of spatial correlations also arises in an- 
other context. Correlation functions are often an important 
statistic for comparing theoretical predictions to the ob- 
serve d distribution of objects. Some an a lytical models exist 
(e.g.. lKaiserlll984l: ICoie fc Kaiserlll989l : iMo fc White] Il996l : 
Shct h fc Tormenll2002l ) that supplement the Press-Schechter 
number densities with additional approximate models, but 
these do not arise naturally within the excursion set ap- 
proach. 

As we review in the following sections, a direct calcu- 
lation of the halo correlation function corresponds math- 
ematically to solving for the simultaneous evolution of two 
correlated random walks with a constant barri e r. Thi s prob- 
lem was first considered by iPorciani et al.l |l998'), who 
made some progress toward a satisfactory solution. We 
l|Scannapieco fc Bark ana 2002) then found an approximate 
but quite accurate analytical solution and used it to calcu- 
late the joint, bivariate mass function of halos forming at two 
redshifts and separated by a fixed comoving distance. We 
showed that our solution leads to a self-consistent expression 
for the nonlinear biasing and correlation function of halos, 
generalizing a number of previous results including those by 
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iKaiseij (|l984h and iMo fe White! (|l996l ). This solution has 
since been used to study, for example, th e impact of clus- 
tered gas minihalos on cosmic reionization (IBarkana fc Loea 
120021 ; lluev et al.ll2005l). inhomogeneous metal enrichme nt at 
high redshift (IScannapieco. Schneider, fc Ferrarall2003l ). and 
observations of metal lin es around Lyman break galaxies 
l|Porciani fc Madaul [20051 ). 

Recently, researchers have found useful applications for 
the more general mathematical problem of ra ndom walks 
with a barrier that is not constant. For instance, ISheth et al.l 
l|200ll ) found that an ellipsoidal collapse model suggests such 
a barrier for defining halos, yielding a model that produces a 
halo mass func tion that better matche s N-body simulations. 
More recently, iFurlanetto et al.l (|2004) used the statistics of 
a random walk with a linear barrier to model the H II bub- 
ble size distribution during the reionization epoch. While 
in principle this distribution could be measured from maps 
of 21-cm emission by neutral hydrogen, upcoming experi- 
ments such as the Mileura Widefield Array and the Low 
Frequency Array are expected to be able to detect ioniza- 
tion fluctuations only statistically, e.g., by measuring the 
correlation function of the 2 1-cm brightness temperatur e 
l|Bowman. Morales, fc Hewittl 1200(1 iMcQuinn et al.ll2006l) . 
While previously approximate expression s for spatial ion- 
ization correlations have b een developed ()Furlanetto et all 
|2004| ; IMcQuinn et~aj] 120051 ). each of these was grafted ex- 
ternally onto the underling formalism, requiring additional 
layers of approximations. 

In this paper w e gen eralize the solution of 
IScannapieco fc Barkana! (|2002h to linear barriers and 
thus develop a self-consistent model for two-point correla- 
tions in this case. The rest of this paper is organized as 
follows. In § 2 we establish our notation and review the 
simplest case of the statistics of a single random walk with 
a constant barrier. We then review in § 3 the generalization 
to a single random walk with a l inear barrier. In § 4 we 
follow the setup and solution of IScannapieco fc Barkana! 
l|2002i ) but generalize it to the case of two correlated random 
walks with linear barriers. Since the barrier corresponding 
to the ionized bubble size distribution during reionization 
is linear to a good approximation, we use this distribution 
in § 5 to illustrate how to apply our results to explore 
various aspects of reionization and of 21-cm fluctuations 
that depend on two-point correlations among the density 
and ionization fields. Our solution, however, is more general 
and can be used in all problems where linear barriers are a 
good approximation to the physical constraints. We briefly 
summarize our results in 5 6. 



2 SINGLE RANDOM WALK WITH A 
CONSTANT BARRIER 

Before considering linear barriers, we first establish our no- 
tation and review in this section the standard derivation of 
the one-point expressions for a constant barrier within the 
contex t of the halo mass function. The basic approach is 
that of lBond et all jl99l|V_ who rederived and extended the 
halo formation model of lPress fc Schechterl l|l974h . 

We work with the linear overdensity field <5(x, 2) = 
p(x, z) I p(z) — 1, where x is a comoving position in space, z 
is the cosmological redshift and p is the mean value of the 



mass density p. In the linear regime, the overdensity grows 
in proportion to the linear growth factor D(z) (defined rel- 
ative to z = 0). The barrier signifies the critical value which 
this linearly-extrapolated 8 must reach in order to achieve 
some physical milestone; in the case of halo formation, an es- 
timate based o n spherical top-hat collapse yields 8 C = 1.686 
(|Peebledll980h in the Einstein-de Sitter model. 

A useful alternative parametrization is to consider the 
linear density field extrapolated to the present time, i.e., 
the initial density field at high redshift extrapolated to the 
present by multiplication by the relative growth factor. In 
this case, the critical threshold for collapse at redshift z be- 
comes redshift dependent even in the Einstein-de Sitter case: 



S c (z) = Sc/D(z) 



(1) 



However, it still represents a constant barrier at any given 
redshift. We adopt this alternative view, and throughout this 
paper the power spectrum P(k) refers to the initial power 
spectrum, linearly-extrapolated to the present (in particular, 
not including non- linear evolution). 

At a given z, we consider the smoothed density in a re- 
gion around a fixed point A in space. We begin by averaging 
over a large mass scale M, or, equivalently, by including only 
small comoving wavenumbers fc. We then lower M until we 
find the highest value for which the averaged overdensity is 
higher than 5 c {z) and assume that the point A belongs to 
a halo with a mass M corresponding to this filter scale. In 
particular, if the initial density field is a Gaussian random 
field and the smoothing is done using sharp fc-space filters, 
then the value of the smoothed 8 undergoes a random walk 
as the cutoff value of k is increased. If the random walk first 
hits the collapse threshold S c (z) at k, then at a redshift 2 
the point A is assumed to belong to a halo with a mass cor- 
responding to this value of k. Instead of using k, we adopt 
the variance as the independent variable: 



Sk= 2^ 



dk' fc' 2 P(fc') 



(2) 



In order to construct the number density of halos in 
this approach, we need to find the probability distribution 
Q(S, Sk), where Q(S, Sk) dS is the probability for a given ran- 
dom walk to be in the interval S to S+d8 at Sk- Alternatively, 
Q(S, Sk) dS can also be viewed as the trajectory density, i.e., 
the fraction of the trajectories that are in the interval 5 to 
5 + d5 at Sk , assuming that we consider a large ensemble of 
random walks all of which begin with 5 — at Sk = 0. 

The distribution Q(S,Sk) satisfies a diffusion equation 



dQ _ 1 d 2 Q 

dSk ~ 2 dS 2 ' 

which is satisfied by the usual Gaussian solution: 



G(S,S k ) 



\Z2ixSk 



exp 



S 2 

■2S k 



(3) 



(4) 



To determine the probability of halo collapse at a red- 
shift 2, we consider random walks with an absorbing barrier 
at S = p, where for halo formation we set v = S c (z). The so- 
lution with the constant barrier in place is given by adding 
an ex tra image solution (|Chandrasekharl 1 19431 : iBond et al.l 

Eaa3): 



Qcon(u, 8, Sk) = G(S, Sk) - G(2v - S, Sk) 



(5) 
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where the subscript "con" refers to the constant barrier 
case. The second ( "image" ) term is clearly (through a simple 
change of variables) itself a solution to the diffusion equa- 
tion, and the combination Q CO n is identically zero on the 
barrier 5 = v, hence it solves the diffusion equation and 
satisfies the required boundary conditions. 

The fraction of all trajectories that have hit the barrier 
v by Sk includes all trajectories except those (represented 
by the solution Q CO n) that still have not been absorbed: 

F>, m n{v, S k ) =1-1 dS Q con (^, 8, S k )=2 / dS G(S, S k 



The differential of this is the first-crossing distribution: 

6 — oo 



r f c \ d P ( q \ f dG(S,S k 

Jcon\t / , Ok) — gg--F>,CDn(,l / , Ok) — I 



3 SINGLE RANDOM WALK WITH A LINEAR 
BARRIER 

The problem of a random walk with a linear barrier has 
been previously considered both in the context of improved 
halo mass functions and the ionized bubble distribution. The 
problem is the same as that considered in the previous sec- 
tion but with a barrier that is linear in the variable Sk, i.e., 
has the form 

S = v + uS k . (13) 

.(6) 

T he first-crossi ng distribution in this case was first derived 
bv lShethl (|l998h . while the full distri bution function Q(6, S k ) 
in this case was first worked out bv lMcQuinn et all f2005). 

A key point that allows our further derivations below 
is that we find a very simple way to express the solution of 
iMcQuinn et~aD (|2005l ): 



2ttS; 



3/2 



exp 



v 

2Sk 



(7) 



where in the second equality we have used the fact that G 
satisfies eq. (|3]). Note that f(v,Sk)dSk is the probability 
that a random trajectory crosses the barrier in the interval 
S k to S k + dS k ■ 

In the halo interpretation, f(v, Sk) dS k is the probabil- 
ity that a given point A is in a halo with mass in the range 
corresponding to S k to S k + dSk- The halo abundance is 
then simply 



dn 
dM 



El 

M 



dS k 



dM 



f{v,Sk) 



(8) 



where dn is the comoving number density of halos with 
masses in the range M to M + dM. The cumulative mass 
fraction in halos above mass M (thus denoted F>) is simi- 
larly determined to be 



F(> M\z) = F > , con (u,S k ) = erfc (n=) 
Note that the complement of this is 



F(< M\z)=F < , can (u,S k )=ert 



nSk 



(9) 



(10) 



While these expressions were derived in reference to 
sharp fc-space smoothing (eq. ([2])), Sk is usually replaced 
in the final results with the variance of the mass M enclosed 
in a spatial sphere of comoving radius r: 



S r (M) = S r (r) 



1 

2^2 



k 2 dkP(k)W 2 (kr) , 



(11) 



where W(x) is the spherical top-hat window function, de- 
fined in Fourier space as 



W(x) 



sin(a;) cos(:r) 



(12) 



The idea of this approach is that the real-space window func- 
tion corresponds more closely to spherical collapse (which 
yielded the critical collapse threshold) , while the mathemat- 
ical problem is simpler in k space and leads to closed-form 
solutions. 



?iin(i / , A*, S,S k ) = G(8, S k ) 



G(2u-5,S k ) 



(14) 



It is easy to check that this simple linear modification of 
the usual image solution of eq. |5]) is identically zero on the 
linear barrier, as required. 

Integrating, we find the fraction of all trajectories that 
have reached the barrier by Sk'- 

d8Q lin (u,a,S,S k ) (15) 

-oo 

V + flSk \ r I V ~ /J-Sk 



erfc 



+ e 



erfc 



<2Sk J \ \f2Sk 

This expression agrees with that in IMcQuinn et all l|2005f ). 
The complement is 

F <t Mniy,ll, Sk) = 1 - F>,lin(^,/i, Sk) (16) 



erfc 



-v - jj.Sk 
V2Sk 



erfc 



v — fJ.Sk 

V2Sk 



Also needed for later is the first moment of the density 
among trajectories that do not hit the barrier: 

5iin(u,fj,,Sk) = / ddSQun(iy,iJ,d,Sk) 



erfc 



>2S k 



(17) 



Finally, we differentiat e to obtain th e first-crossing distribu- 
tion in agreement with IShethl (|l998l ): 

d 

flinty, Sk) = -^-F >t Un(.V,fJ>,Sk) 



2-kS 



3/2 



exp 



(v + fiSkf 



2Sk 



(18) 



4 TWO CORRELATED RANDOM WALKS 
WITH LINEAR BARRIERS 

4.1 Analytic Preliminaries 

We follow Scannapi eco fc Barkanal (|2002l ) in setting up the 
problem of the statistics of two correlated random walks. 
We consider points A and B separated by a fixed comov- 
ing distance d. Note that this definition of distance is in 
Lagrangian space, which is intrinsic to any Press-Schechter 
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type approach. Thus, it is the comoving distance between 
points A and B at early times, and does not take into ac- 
count subsequent peculiar motions of these points. If we con- 
sider smoothed densities identified by sharp fc-space filters 
fci at point A and fc 2 at point B, then the cross-correlation 
of the densities involves only those k values common to both 
filters, and its value is 



k , 2dk , S m(k'd) , 
k'd 



(19) 



where the upper integration limit is k = min[fci,fc 2 ], and 
when we write £k as a function of Sk it is related to k by 
eq. (J3]l. It is also convenient to define 



so that 



sin [fc(S fc ) dj 
k(Sk) d 



V (d,S' k )dS' k 



(20) 



(21) 



Just as the real-space variance is often used in the one- 
point case, the two-point quantities we discuss below will 
use the correlation between two spatial filters centered about 
two points at a separation d. In this case, the standard ex- 
pression is 



e r (d,n,ra) = ^ 



k 2 dk sm( ft rf ) p(k)W (kn)W (kr 2 ) 
kd 



(22) 



where ri and r 2 are the radii of the two filters, and W(x) is 
again the top-hat window fun c tion g iven by eq. (fT2")l . How- 
ever, IScannapieco fc Barkanal (|2002T ) showed that in order 
to ensure the most physically-reasonable behavior of the so- 
lution in various regions of the parameter space (including 
in limits that reduce to the one-point case), it is better to 
substitute for £ k the quantity 



ir maa {d, n,r 2 ) = £ r [d, max(ri, r 2 ), max(n, r 2 )] , 



(23) 



which is equal to £ r when the two filters have equal radii. 
When the variances Si and S2 are used as the fundamental 
variables, we find the corresponding r\ and r 2 by inverting 
the relation S r (r) given by eq. (|11[) . 



4.2 Basic Setup 

We continue to follow IScannapieco fc Barkanal (I2002T ) as 
we consider simultaneous correlated random walks of two 
overdensities 8i(S k ,i) and <5 2 (Sfc i2 ) separated by a fixed La- 
grangian distance d. As in the one-point case, for the deriva- 
tion we adopt sharp fc-space filters. We want to determine 
the joint probability distribution of these two densities, 
Q(6i, 82, Sk,i, S k ,2,d). In terms of a trajectory density in the 
(Si, S2) plane, Q(5i, 82, Sk,i, Sk,2, d) d8i d<5 2 is the fraction of 
trajectories that are in the interval 81 to Si + d8i and 82 to 
5 2 +d<5 2 at (Sk,i, Sk,2)- Below we will take Sk,i and Sk,2 to be 
the final variances of these trajectories, denoting interme- 
diate variances with the primed notation S' k 1 and S' k 2 . We 
then consider a large number of random walks all of which 
begin with 81 = and 82 = at S' kl = and S' k 2 = 0. 

With sharp fc-space filters, the problem simplifies due to 
the fact that we are working with a Gaussian random field. 



IScannapieco fc Barkanal (|2002T ) showed that we can consider 
Q to be a function of a single variable S' k , with a diffusion 
equation 



( 1 g2 g I n(d g"| B 2 Q . 1B 2 Q 
2 ~ ''V"' a k) dS 1 S 2 ~ 2 



8Q 



= < 



2 86^ 

1 d 2 Q 

2 dS'f 

1 d 2 Q 

2 W 



^k ^ ^ k ,min 
Sk,2 < S'k < Sk,l 
Sk.l < S'k < Sk,2 

(24) 



where Sk mm is the smaller of S k 1 and Sk 



4.3 Two-Step Approximation 

While the full solution of the double barrier problem r equires 
a numerical approach. IScannapieco fc Barkarial (1200 21 ) found 
a simple approximate analytic solution that captures the 
underlying physics of two-point collapse. 

Consider the expression for the differential correlation 
coefficient r/(d, S' k ), eq. (|20[) . While this is an oscillating func- 
tion, it equals unity at small values of S' k and its amplitude 
declines towards zero once kd ^> 1. Thus, for small S' k val- 
ues, the two random walks are essentially identical, while at 
large S' k , the two random w alks become independent . 

These observations led fScannapi eco fc Barkanal (|2002T l 
to propose a "two-step" approximation in which ri(d,S' k ) is 
replaced with a simple step function. In order to preserve 
the exact solution for Q at S' k = Sfc, m m in the absence of the 
barriers, we specifically took 



v(d,S' k ) 



1 0< SK&(d,Sfe, min ) 



£ h (d,S kll 



<S' k ^ S k 



(25) 



Hereafter we adopt a general notation for the variances 
and correlation functions, using S to represent either the 
fc-space filtered quantity, Sk, its real space equivalent S r , 
or any alternative definition. Similarly, £ denotes £r max or 
^k(r)- Following the common approximation taken in the 
single-particle case, in all applications we use the real-space 
quantities, i.e., Si and 5 2 denote S r (Mi) and 5V(M 2 ), re- 
spectively. Also, we use for £ the correlation function £r max 
as given by eq. 1|23[1 . Note that although we write the de- 
pendence of various functions on £ explicitly, £ is not an 
independent variable but instead is a function of Si and S2 
(as well as the separation d). 



4.4 Analytic Solution With Linear Barriers 

Using the two-step approximation, IScannapieco fc Barkanal 
(|2002r ) found the analytic solution for Q in the case of con- 
stant absorbing barriers at 81 = vi and 82 = vi- We follow 
their derivation, but generalize it to the case of two (possibly 
different) linear barriers at 81 — vi+fiiS' and 82 = V2+H2S' . 
Under the two-step approximation, we must first evolve Si 
for 0^5*'^^. Since we are assuming that the two random 
walks are identical in this regime, we must place the barrier 
on Si at Si = v m + HmS' where we choose ^ m and p m to be 
as large as possible such that the resulting barrier still lies 
below both of the original linear barriers, throughout the 
relevant range of 5*'. In principle, the best approximation 
would be to adopt at each S' the lower of the two barri- 
ers. However, if the barriers were to cross within the range 
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^ S' ^ £, this would require an extra convolution com- 
pared to our solution below and would thus complicate it 
substantially. Fortunately, this appears not to be needed in 
practice, at least in our main application which is reion- 
ization. In the examples given in § 5, we find that while v 
changes rapidly with redshift, /i (which does change in the 
opposite direction) varies extremely slowly, so that if barriers 
are considered at two different redshifts, the lower-redshift 
one is the lower barrier at all relevant values of S' . 

Quantitatively, the solution for a single linear absorbing 
barrier, eq. (I14|l . gives Q at 5' — £: 



Qa(Vm, Jim, <5l , S 2 , £) = 



-2h> m fi n 



G(2v m -8i,£)] 



x 5 D {Si - S 2 )0(v m + Mm£ - 81) 



(26) 



where So is a one-dimensional Dirac delta function and 9 
is the Heaviside step function. We then evolve the random 
walks in 81 and 82 independently from their common start- 
ing point at £ up to 5i and S2, with the barriers at vi + Mi5' 
and V2 + H2S' , respectively. Thus, we first convolve eq. (f26|) 
with the no-barrier solutions for the two independent ran- 
dom walks, 

Q b (8i,8 2 ,Si,S2,0 =G(8i,Si -S)G(6 2 ,S 2 - £) • 

(27) 

f , we can write this 



Letting 8 be the value of 81 at S' 
convolution explicitly as 

Qo(vm, Mm, 81, 82, Si, S2, £) = 

Vm+Puit 

d8 Q\i u {Vm, fha, 8, £) 



x G(8i -8,Si -£)G(6 2 -8,S 2 - f ) ■ 
Evaluating this yields 

Qo(Vm, Mm, Si, 82, Si, S2, £) = 

Q+{Vm + Mm£, 81,82, Si, 52, C) 

+ e 

x Q-(v m — /i m ^,2^ m — 8i,2v m — 52,5*1,52,0 

where 

Q±(v, 81,82, Si, S 2 ,0 
1 



(28) 



(29) 



4^^5152 -e 

8iS 2 + 82S1 — 28182^, 



x exp — 




(30) 



and we have defined 



s 



5i5 2 -C 2 

v 81 82 

S ~ 5i - l~ 5 2 - £ 



(31) 



Note that the quant i ties Q± are unchanged from 
IScannapieco k, Barkanal l|2002r i. but the solution for Qo is 
now more general. 

Finally, we must account for the additional barriers 



on Si and 82 in the regime where their random walks are 
independent. To do this, we first note that the barrier 
81 = vi + niS' can be written also as a linear barrier in 
terms of the relative variables Si — S and 5' — £: the bar- 
rier is at Si — 8 = [vi + mi£ — 8] + /ii(5' — £). Thus, the 
linear-barrier solution of eq. l|14p shows that we must sub- 
tract from the no-barrier term G(Si — 8, Si — £) in eq. 
an image-like second term: 



-2(fi+Ml£-6)M 



1 G{2{vi + mi£ - 8) - (Si - 8), Si-0 ■ (32) 



Thus, the solution Q can be written as 
Q(vi,v 2 ,^i,^ 2 ,8i,8 2 ,Si,S 2 ,S) = 

<Aii + MmC 

d8 Qlin(l^m,Mm,(5,0 



x Q\\n{vi + fJ.i£ — 8, mi, 81 — 8, Si — 

x Qun^ +M -8,^2,82- 8, S2-O ■ (33) 

This integral is difficult since the exponential factor 
that multiplies G in eq. (|32p itself contains the integration 
variable 8, and the result therefore cannot immediately be 
written in terms of Qo- However, we solve this difficulty by 
noting that the expression in eq. (|32|> can be written in an 
equivalent, alternate form: 

e -2(5 1 -(„ 1+M1 S 1 ))Ml G((5i +S _ 2{yi + Si _ i (34) 

Thus, the solution Q can also be written as 
Q(vi,v 2 ,^i,^ 2 ,8i,8 2 ,Si,S 2 ,S,) = 

<Aii + MmC 

dS Qlin(l / 111 • Mm, <5,£) 

x Qun(5i - (ui + mi5i),mi,5i -6, Si - £) 

x Qh„(5 2 - (v2 + tM 2 S2),H2,82-8,S 2 -0 . (35) 

This integration yields the complete solution: 
Q(vi,v 2 , Mi, M2, 81,82, Si, 52,?) = 

Qo(^m, Mm, Si, 8 2 , Si, S2, C) 

+ exp [2(8^ - 8i)fii + 2(8? - 82)^2} 
x Qo(fm,Mm,25i r — 8i,28\ v — 8 2 ,Si,S 2 ,S) 

- exp [2(S 2 r - S 2 )n 2 ] 

x Qo(u m , ^i m , Si, 2S 2 ' — 82, Si, S 2 , 

- exp [2(<5i r - <5i)mi] 

x Qo(^m, Mm, 2<5i r — 8i,S 2 ,Si,S 2 ,£,) , (36) 
where we have defined the 8 values on the barriers: 



Si 1 = vi + mi5i; 52 r = ^2 + M252 



(37) 



IScannapieco fc Barkanal (|2002l ') showed that the solu- 
tion with the two-step approximation is very accurate in the 
case of constant barriers, giving the same results as a full nu- 
merical solution to within at most 2%, with the difference 
typically much smaller than this value. Since the idea of the 
approximation (as presented in the previous subsection) is 
based on the properties of two correlated walks and not on 
any particular property of the barriers, we expect this ap- 
proximation to be accurate in the case of linear barriers as 
well. 
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4.5 Bivariate Cumulative Distribution 

Having developed in the previous subsection an accurate 
approximation to the joint statistics of two correlated ran- 
dom walks, we now apply this distribution to find the joint 
probability of having the two random walks cross their re- 
spective barriers before reaching two given values Si and 5*2. 
In particular applications, this quantity can be interpreted 
as the joint probability that point A is in a halo above a 
mass Mi(Si) and point B is in a halo above a mass A^S^), 
or as the joint probability that point A is in an ionized bub- 
ble above some size (see § 5.1) and point B is in an ionized 
bubble above some other given size. 

Consider first the following quantity: 

F < (vi,V2,(Ll,(L2,Sl,S2,V) = 

dSi I dS 2 Q(yi,V2, m, (12,61,82, Si, S2,£) 



(38) 

This is the probability that both random walks are not ab- 
sorbed before reaching the point (Si, 82)- We denote it F< 
since, e.g., in the halo-formation case it is the chance that 
point A is in a halo below a mass Mi (Si) and point B is in 
a halo below a mass M2(S2), 

We can find a n expression for F < in t e rms o f a sin- 
gle integral (as did Scannap ieco fc Barkanal (|2002T ) in the 
constant-barrier case), by writing Q in the form of eq. (|33[) 
and performing the Si and 82 integrals. The result is 

F < {vi,v 2 , fii, (J.2, Si,S 2 ,0 = 

dS Q\i n (v m , jU m , <5,£) 



X -F<J in (^l + (ll£ — 6,(1,1, Si — £) 
X F < s ln (u2 + (12£, — 5,(I2,S2 - f) 



(39) 



This expression is easy to understand: After the correlated 
random walk reaches 8 at S' — f, without hitting the joint 
barrier v m + /U m S', the subsequent random walks are inde- 
pendent. We therefore multiply the probability that random 
walk #1 does not hit its barrier between S' = £ and Si, with 
the probability that random walk #2 does not hit its bar- 
rier between S' = £ and S2. This is then integrated over 
the probability distribution of reaching various values of 8 
at S' = £. 

The complementary quantity f> is the probability that 
both random walks are absorbed before reaching the point 
(Si,S2). This cannot be calculated with a similar expres- 
sion as in eq. just replacing F <i \ in with F >t u n in the 
integrand, since the barrier can also be crossed before the 
correlated walk reaches S' — £. Instead, we find f> as the 
complement of the chance that at least one of the random 
walks is not absorbed. The latter chance equals the chance 
that #1 is not absorbed, plus the chance that #2 is not 
absorbed, minus (to eliminate double counting) the chance 
that both are not absorbed. We obtain from this: 

F > (vi,V2,(J.l,(l2,Sl,S2,£,) = 

1 + F < (lS 1 , V 2 , (ll, (12, Si, S2,$.) 

-F < ,lin(»l,(J-l,Sl)-F < ,l ill (v2,IJ l 2,S2) . (40) 

We can similarly calculate the mixed quantities; e.g., the 



chance that walk #1 is not absorbed but #2 is absorbed is 

F <> {V1,V2, (J-l, (12, Si, S 2 , £) = 

i ? >,lin(^2,M2,S 2 ) - F > (ui, V2,(ll,(J,2,Sl,S2,i) ■ (41) 

Finally, the chance that walk #1 is absorbed but #2 is not 
absorbed is obtained by switching the indices 1 and 2 in 
eq. (|HJ. 

4.6 Other Distributions 

The bivariate first-crossing distribution / dSi dS2 is the 
probability of having random walk #1 cross the barrier in 
the range Si to Si + dSi and point #2 cross in the range S2 
to S2 + dS2. This is simply related to the bivariate cumula- 
tive distribution as 



f(Vl,V2,(ll,(l2,Sl,S2,0 
!) Q 

-F < (vi,V2, (11,(12, Si, S2, 



(42) 



<9Si dS 2 

where £ is not considered an independent variable (and 
so the partial derivatives involve variations of £). The ex- 
pression for / can in principle be simplified by bringing 
the derivatives inside the integral in eq. (|39p and using 
the properties of the int e grand as in the analogous case in 
IScannapieco &: Barkanal (|2002l ). although here the expres- 
sions are more complicated (and there is also a contribution 
from the £ that appears in the integration limit). Since we 
are not directly interested in / in the context of upcoming 
probes of reionization, we do not develop this further here. 
The derivatives in eq. (|42[l can also be evaluated numerically. 

Various correlated distributions of density and of ion- 
ization (i.e., hitting the barrier) can also be calculated with 
our solution. For example, consider the probability distribu- 
tion of Si at Si given that random walk #1 has not been 
absorbed by its barrier while #2 has been absorbed by its 
barrier before S2. We calculate this as follows: After the 
correlated random walk reaches 5 at S' = £, without hitting 
the joint barrier f m + ^ m S' (so that #1 will be unabsorbed), 
the subsequent random walks are independent. We therefore 
multiply the probability that random walk #1 reaches Si at 
Si without hitting its barrier on the way, with the prob- 
ability that random walk #2 does hit its barrier between 
S' = £ and S2. This is then integrated over the probability 
distribution of reaching various values of S at S' = £. The 
probability is proportional to the following quantity: 

f[5 1<>] (^l) v i, Mi> ^2, 61 , Si, S2, = 

dS Ql in (Vm, (ljn,S,$,) 



x Qiin(<5i — (vi + (uSi), (ii,Si — S, Si - 

x F > , Un {v2 + (i2i - 5, (12, S2-O ■ (43) 

As written, this distribution for 81 is not normalized; a nor- 
malized probability distribution can be obtained by dividing 
by the probability in eq. (|4ip (with indices switched) that 
walk #1 is not absorbed by its barrier while #2 is absorbed. 

Finally, consider the probability distribution of <5i at 
Si given that both walks have not been absorbed. This is 
similar to eq. (|39|l except that we integrate only over the 
values of #2. Thus, we first calculate the quantity: 

f[6\<](vi,V2,(tl,(l2,5l,Sl,S2,0 = 
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/ dS Qlin(Vm, A*m, S, f) 

J 5— — 00 

x Qnn(8i — {vi + mSi),m,5i — S, Si — £) 

x F <tlia (y2+H2£-6,H2,S2-ti) ■ (44) 

A normalized probability distribution can be obtained from 
this by dividing by the probability in eq. (|39|l that both 
walks are not absorbed. 



5 ILLUSTRATION: COSMIC REIONIZATION 

5.1 The Density and Ionization Fields 

We illustrate the a pplication of our solutio n to reionization 
using the model of iFurlanetto et"all (|2004l ') for the ionized 
bubble distribution. According to this model, a given point 
A is contained within a bubble of size given by the largest 
surrounding spherical region that contains enough ionizing 
sources to fully reionize itself. If we ignore recombinations, 
then the ionized fraction in a region is given by C/coii, where 
/coil is the collapse fraction (i.e., the gas fraction in galac- 
tic halos) and £ is the overall efficiency factor, which is the 
number of ionizing photons that escape from galactic halos 
per hydrogen atom (or ion) contained in these halos. This 
simple version of the model remains valid even with recombi- 
nations if the number of recombinations per hydrogen atom 
in the IGM is treated as uniform; in this case, the resulting 
reduction of the ionized fraction by a constant factor can be 
incorporated into the value of £. 

In the extended Press-Schechter model [compare 
eq. @] , in a region containing a mass corresponding to vari- 
ance Sm> 



- erfc 



\ \J 2(tS'miii Sm ) 



(45) 



where S m in is the variance corresponding to the minimum 
mass of a halo that hosts a galaxy, and 8 m is the mean 
density fluctuation in the given region. While this describes 
fluctuations in / co ii well, the cosmic mean collapse frac- 
tion (and thus the overall evolution of reionization with 
redshift) is better desc ribed by the halo mass function of 
ISheth fc Tormenl ll 19991) (wi t h the updated parameters sug- 
gested by ISheth fc Tormenl (|2002h ). We thus use the lat- 
ter mean mass function and adjust / co ii in different re- 
gions in proportion to th e extended Press-Schechter formula; 
iBarkana fc Loebl l|2004l ) suggested this hybrid prescription 
and showed that it fits a broad range of simulation results. 
The resulting condition for having an ionized bubble of a 
given size, writ ten as a condition for S m vs. S m , is of the 
same form as in IFurlanetto et all ([2004), at a given redshift, 
and thus (as they show ed) yields a linear barri er to a good 
approximation (see also Furlanctto ct al. (200ji|)). 

In the model of IFurlanetto et aL I (|2004n T the total frac- 
tion of points contained within bubbles, as given by the 
model [i.e., eq. (|15p ], comes out slightly different from the 
direct result for the mean global ionized fraction, Xi = C/coii 
in terms of the cosmic mean collapse fraction. To deal with 
this, we adopt the direct values of Xi versus redshift (or the 
values measured in a simulation, when comparing to one), 
and adjust £ within the model to an effective value of C, at 



x 

A 




d [com Mpc] 

Figure 1. Ionization correlations. We show the joint probabil- 
ity -F> that two points separated by a distance d are in ionized 
regions, divided by the mean ionized fraction We consider 
Xi = 0.1, 0.2, 0.4, 0.6 and 0.8 (bottom to top), assuming param- 
eters for which the universe fully reionizes at z = 6.5 due to stars 
in halos with efficient atomic cooling. 



each redshift that gives a model value of Xi that equals the 
desired one. 

We illustrate the power of our solution from § 4 by cal- 
culating a number of different statistics. In the f ollowing ex- 
ample s we use the cosmological parameters from lZahn et al.l 
(2006) since we compare with their results in the following 
subsection. In this subsection we assume that the efficiency 
£ is constant in all halos with circular velocity V c above 
16.5 km/s (corresponding to efficient atomic cooling); let- 
ting reionization end, e.g., at z = 6.5, yields a real £ = 8.9 
(which is held fixed, independent of redshift, unlike the ef- 
fective model Q. Figure 1 shows the probability that two 
points are both in ionized regions, divided (for visual clar- 
ity) by the mean ij, asa function of the distance between the 
points. The probability is F> as given by eq. (|40[) evaluated 
at Si = & = Smin- This probability in our model naturally 
satisfies the limits F> — > Xi when d — > (perfect correlation) 
and J> — > xf when d — > 00 (no correlation) , while these lim- 
its had to be artificially inserted into previous models for the 
ionization correlation function. The characteristic distance 
at which F > makes the transition between these two lim- 
its grows as reionization proceeds, reflecting the increase in 
the characteristic bubble size as larger and larger groups of 
galaxies produce joint ionized regions. 

Our solution also allows us to calculate density- 
ionization correlations. Figure 2 shows the probability distri- 
bution of the density fluctuation on the scale S m i n around a 
point A. In this subsection and the next, we use the notation 



5i(zi) = D(zi)5i 



(46) 



where the growth factor converts from the linearly- 
extrapolated Si at redshift (which we have been using) to 
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Figure 2. Density-ionization correlations. We show in each panel 
the normalized probability distribution of the density fluctuation 
on the scale of 0.09 com Mpc (corresponding to S mln = 46) 
around a point A, in general (dotted curve) or given that point 
A is in a neutral region (solid curve). For the latter, we also show 
the break-down into a contribution from the case where point B 
is ionized (short-dashed curve) or where point B is neutral (long- 
dashed curve). With the same assumptions as in Figure 1, we 
consider X{ = 0.8 (which implies z = 7.3) and a separation d = 1 
com Mpc (top panel) or 10 com Mpc (bottom panel). 



the linearly-extrapolated 8i at redshift z\ . Given that point 
A is neutral, we calculate as detailed in § 4.6 the separate 
probability distributions of 5i{zi) given that a point B a 
distance d away is either ionized or neutral. When d = 1 
com Mpc the two points are highly correlated, and point B 
is most likely to be neutral as well, especially when 5i(zi) 
is very negative. However, when d = 10 com Mpc the corre- 
lations are weaker, and point B is most likely ionized, but 
only with a 65% chance although the IGM as a whole is 80% 
ionized in the plotted example. 



5.2 The 21-cm Power Spectrum 

During cosmic reionization, we assume that there are suffi- 
cient radiation backgrounds of X-rays and of Lya photons 
so that the cosmic gas has been heated to well above the 
cosmic microwave background (CMB) temperature and the 
21-cm level occupations have come into equilibrium with the 
gas temperature. In this case, the observed 21-cm brightness 
temperature relative to the CMB is independent of the spin 
temperatur e and, for our assum ed cosmological parameters, 
is given by l|Madau et al.lll997l ) 



T b = T 6 *; T b = 25 



1+z 



8 



mK 



(47) 



with ^ = thi[1 + S(z)], where xni is the neutral hydrogen 
fraction and we also used the notation of eq. (|46[) . Under 
these conditions, the 21-cm fluctuations are thus determined 



by fluctuations in ty. To determine its statistical properties 
using our model, we note that for a given random walk, the 
value of the neutral hydrogen is either 1 (if the barrier has 
not been pierced) or (if it has). 

Consider points A and B at a distance d from each 
other. Then the correlation function of <!< is = (^1^2) — 
(5 , i)($2), where the mean value is, e.g., for point 1: 



= F <> u n (l/i, (Ml, Si) + D(zi) 5u n (fl,/Lti, Si) 



(48) 



where we set 5*i = S m i n and use eq. (|17[1 . This result arises 
from the fact that the average value of 9 is simple the value 
of [1 + S(z)] averaged only within neutral regions (where 
xm = 1, corresponding to random walks that have not been 
absorbed); e.g., the first term (unity) yields simply the frac- 
tion of the universe which is still neutral. 

To average over the value of [1 + 5i(zi)] x [1 + ^2(22)] 
when both points are neutral, we follow the derivation of 
eq. (OH, writing, e.g., (1 + D(zi)5i) = (1 + D(zi)S) + 
D(zi) (5i — S). Note that our solution for Q describes ex- 
actly those random walks that correspond to both points 
being neutral (i.e., not absorbed by the barrier). We obtain 



d5 Qn n {v .1:1:1. « 

— CO 

(1 + D{zi)5)F <Ml ^i + fj.it - S,m,Si - 
+ D(zi) <5iin(fi + a*i£ — 8,m,Si - £) 
(1 + D(z 2 ) 5)F <Xin (v 2 + a*2? - 5, /j,2, S 2 



+ D(z 2 ) <5lin(f2 + — <5, M2, 5*2 — £) 



■0 

(49) 



We use this equation with Si = S 2 — Smm. However, 
we can calculate the correlation function of ^ down to 
smaller scales by assuming that the smaller-scale power 
does not influence the ionization and is independent of the 
larger scales. Thus, we simply add to the contribution 
D(zi)D(z 2 )[£ r (d,0,0) - £ r (d, r min ,r min )] times the proba- 
bility that both points are neutral, where r m i n is the scale 
corresponding to Smm- Note, though, that these scales are 
quite small and we expect large non-linear corrections to the 
model over this range of scales. 

Having calculated the correlation function of the 21-cm 
brightness temperature, we Fourier transform it and obtain 
the power spectrum P(k) as a function of wavenumber k. 
We express the result in terms of a characteristic quantity 
that has units of temperature: 



A 2 i(fc) = f b v /PP(fc)/(2 7 r 2 ) 



(50) 



Figure 3 shows that the 21-cm power spectrum changes 
shape during reionization, acquiring large-scale power and 
flattening on scales up to several tens of Mpc as the char- 
acteristic bubble size grows towards the end of reionization. 
The amount of large-scale power depends strongly on the 
bias of the typical ionizing galaxies. The bias is larger for 
more massive halos, leading to stronger large-scale fluctua- 
tions in this case. These trends are qualitatively similar to 
those seen in previous approx imate models that were con- 
structed more artificially [e.g., iFurlanetto et al.l (|2006f )]. 

As a final example, in Figure 4 we compare our model 
quantitatively to a numerical N-body plus radiative trans- 
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Figure 3. 21-cm power spectrum. Assuming parameters for 
which the universe fully reionizes at z = 6.5, we show in each 
panel results for Xi = 0.1 (dot-dashed curve), 0.3 (dotted curve), 
0.5 (long-dashed curve), 0.7 (short-dashed curve) and 0.9 (solid 
curve). We consider stars forming in all halos above V c = 16.5 
km/s (top panel; corresponds to efficient atomic cooling) or only 
in halos ten times as massive, above V c = 35.5 km/s (bottom 
panel; corresponds to strong feedback in low-mass halos, e.g., due 
to photoheating or supernovae). 



fer simulation bv lZahn et all l|2006l ). In this comparison we 
modify the model slightly in acco rdance with the a ssump - 
tions in their simulation. Unlike Furlanetto et all l|2004l) . 
who effectively assumed that the star formation rate in halos 
is prop ortional to the rate of gas infall into them. lZahn et all 
(2006) assumed a constant mass-to-light ratio, which sets 
the star formation rate in halos to be proportional to their 
total gas content at a given time. This assumption leads 
to a slightly different condit ion for havi n g eno ugh sources 
to ionize a given region [see IZahn et al.l (2006)], which we 
again approximate as a linear barrier constraint. We set the 
minimum halo mass to be 2 x 1O 9 M0, as assumed in the 
simulation, and set the effective efficiency factor at each 
redshift so that the model yields the same global ionized 
fraction as measured in the simulation. We also compare 
our r esults to the nume rical extended Press-Schechter model 
from lZahn et al.l l|2006l ). where they numerically applied the 
spherical ionization condition (in real space) to the linear 
density field. 

These comparisons are an ambitious challenge for our 
model since it is fully analytical and makes necessary ap- 
proximations in using spherical averages in the statistics, in 
applying simplifying assumptions that are strictly valid only 
in fc-space, and in neglecting significant non-linear correc- 
tions. The model also relies on the two-step approximation, 
approximates the reionization condition as a linear barrier, 
and is based on a Lagrangian approach. The simulation is 
limited as well, with fluctuations in the measured A2i(fc) 
indicating a lack of convergence on large scales, while on 



Figure 4. 21-cm power spectrum. We compare our model predic- 
tion (solid curves) to those from the simulation (dashed curves) 
and numer i cal ex tended Press-Schechter (dotted curves) from 
IZahn et al.l J2006h . The results are shown at several different red- 
shifts, as indicated in each panel. At each redshift z we adjust the 
value of the efficiency in our model in order to match the mean 
global ionized fraction Xi from the simulation. 

small scales the mass resolution corresponds to only 64 par- 
ticles per 2 x 10 9 Mq halo, well below the 500 required for 
reaso nable confidence as indicate d by careful convergence 
tests (S pringel fc Hernq uist 2003). Nevertheless, the com- 
parison indicates that the simple analytical model captures 
the correct trends such as the change in power-spectrum 
shape with redshift, and can therefore be used to estimate 
the quantitative results and to explore the dependence on 
model parameters such as the astrophysical properties of the 
ionizing sources. 



6 SUMMARY 

We have presented an approximate but fairly accurate an- 
alytical solution to the mathematical problem of the joint 
evolution of two correlated random walks with linear absorb- 
ing barriers. Our self-consistent sol ution is a generalization 
of the constant-barrier solution of IScannapieco fc Barkanal 
(2002) and is based on their two-step approximation. Physi- 
cally this mathematical setup can be applied to a number of 
topics in galaxy formation where spatially-dependent feed- 
back or two-point correlations are important. 

We have emphasized in particular the direct relevance 
to extended Press-Schechter models of t he ionizing bub- 
ble d istribution during cosmic reionization (jFurlanetto et all 
|2004| ). In this context, the joint probability distribution Q 
of the random-walk trajectories [eq. (|36[) ] corresponds to the 
bivariate density distribution at two points when both points 
are neutral. The bivariate cumulative probability F > of both 
points hitting their barriers [eq. (|40|l ] corresponds to the 
probability that both points are in ionized regions. Other 



10 R. Barkana 



distributions [e.g., as given by eq. (|43j] correspond to var- 
ious elements of the joint correlations among the densities 
and ionization states of the two points. 

We have shown that our model can be used not only to 
calculate density-ionization correlations (Figures 1 and 2), 
but also (Figure 3) the power spectrum of fluctuations in the 
21-cm temperature brightness, which may be observed in the 
next few years. Like any analytical approach to complicated 
non-linear physics, our model is approximate and simplified 
in a number of ways, but it correctly captures the trends 
seen in simulations of reionization (Figure 4) and thus can 
be used to explore various scenarios of cosmic reionization 
and their observable consequences. 
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